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Abstract. It is known that in gravitational instability 
ON scenarios the nonlinear dynamics induces non-Gaussian 
features in cosmological density fields that can be in- 
P vestigated with perturbation theory. Here, I derive the 
expression of the joint moments of cosmological density 
fields taken at two different locations. The results are 
lf~) valid when the density fields are filtered with a top- 
hat filter window function, and when the distance be- 
tween the two cells is large compared to the smoothing 
length. In particular I show that it is possible to get 
the generating function of the coefficients C VA defined by 
^ (8P(x 1 )8"(x 2 )) c = C Ptq (8 2 (x)) p+q - 2 (8(x 1 )8(x 2 )) where 
f^) 6(x) is the local smoothed density field. It is then possi- 
^sO ble to reconstruct the joint density probability distribution 
function (PDF), generalizing for two points what has been 
Q obtained previously for the one-point density PDF. I dis- 
Q^cuss the validity of the large separation approximation in 
q an explicit numerical Monte Carlo integration of the 6*2,1 

parameter as a function of \xi — x 2 \- 
CZJ A straightforward application is the calculation of the 
^ large-scale "bias" properties of the over-dense (or under- 
dense) regions. The properties and the shape of the bias 
function are presented in details and successfully com- 
pared with numerical results obtained in an N-body sim- 
ulation with CDM initial conditions. 

Key words: Cosmology: theory - large-scale structure of 
the universe - Galaxies: clustering 



1. Introduction 

Recently a lot of progress has been made in the study of 
the nonlinear features induced by the gravitational dy- 
namics in an expanding Universe. In particular it has 
been shown that it is possible to apply the techniques 
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of perturbation theory (PT) to derive statistical proper- 
ties of the density field at large-scale. For instance it is 
possible to use PT to calculate the leading order of the 
low order reduced moments of the local density, (Pee- 
bles 1980, Fry 1984). Refining these calculations Goroff 
et al. (1986) showed that the unavoidable filtering of the 
cosmological density fields affects the results. They calcu- 
lated numerically those effects in a particular case. More 
recently Juszkiewicz (1993), Bernardeau (1994b), Lokas 
et al. (1995) investigated by means of analytic calcula- 
tions these smoothing effects either for a top-hat or a 
Gaussian window function. And finally the initial results 
obtained for the third and fourth cumulants of the one- 
point density distribution function, have been extended to 
the whole series of the cumulants by Bernardeau (1992) 
when the smoothing effects are neglected and Bernardeau 
(1994c) for top-hat filtering. These results are found to be 
in remarkable agreement with results of numerical simula- 
tions (Juszkiewicz et al. 1995, Bernardeau 1994b, c, Baugh, 
Gaztanaga & Efstathiou 1995, Gaztanaga & Baugh 1995, 
Lokas et al. 1995). So far, however, the properties of the 
moments of the two-point density PDF have never been 
investigated. 

This is however a quantity of various interests. Its most 
direct application would probably be for a proper estima- 
tion of the errors due to finite sampling for the measure- 
ments of the statistical properties of the one-point den- 
sity field. Indeed, when the density is measured in a finite 
sample at different locations, these measures are not inde- 
pendent and the estimation of the errors should take into 
account these correlations. The two point density PDF is 
expected to provide the most dominant part of these cor- 
relations (Colombi, Bouchet & Schaeffer 1995, Szapudi & 
Colombi 1995). 

The derivation of the joint-density PDF is also of cru- 
cial interest for the study of the properties of a gravi- 
tational fluid. In particular it can give insights on how 
the dense spots of the field are correlated together, com- 
pared to the global matter correlation function. This is 
thus related to the general problem of the bias in large- 
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scale structure formation. The problem that shall even- 
tually be addressed here is the estimation of the gravity 
induced bias in the quasi-linear regime. 

The most commonly considered bias factor is the bias 
b associated to the galaxies. In general the problem of the 
bias emerges because the mass distribution of the universe 
cannot be observed directly and thus has to be estimated 
with the apparent distribution of tracers. However, these 
tracers, galaxies, (or galaxy clusters, super-clusters) may 
not be reliable tracers of the underlying mass distribution 
so that their density fluctuations do not reproduce the 
underlying mass density fluctuations. From the observa- 
tions of the density fluctuations of the clusters and the 
galaxies, it is already clear that there is a relative bias be- 
tween these objects. This concern led people to introduce 
a bias parameter b relying the distribution of tracers to 
the mass distribution. In its simplest form it is assumed 
that the local (smoothed) over-density of objects, 6 ga i., is 
proportional to the local mass over-density, 

<5gai. = b 6 mass . (1) 

The physical processes that can possibly affect the value of 
b are numerous, but at very large scale it is likely that only 
gravity can substantially affect the correlation of large 
mass concentrations. Neglecting, at first sight, the prob- 
lem of identifying the actual astrophysical objects in the 
evolved final density field, we are led to a purely grav- 
itational physics problem, that is the calculation of the 
correlation of dense spots in a self-gravitating fluid. 

There have been attempts to solve this latter prob- 
lem. For instance Doroshkevich & Shandarin (1978a, b) 
and Bardeen et al. (1986) derived the bias property of the 
peaks in a Gaussian density field. This idea had been sug- 
gested by Kaiser (1984) to explain the relatively strong 
clustering of the Abell clusters. This gave a good insight 
of the possible selection effects in a random field. This 
result, however, did not take into account the subsequent 
nonlinear dynamics which is likely to affect the correlation 
properties of the peaks, such as their peculiar displace- 
ment, or more worrying, the possible merging of halos. 
In an attempt to overcome these difficulties and start- 
ing with a different point of view - a phenomenological 
description of the matter p-point correlation functions in 
the final nonlinear field - Schaeffer (1985) and more re- 
cently Bernardeau & Schaeffer (1992), extending a work 
of Balian & Schaeffer (1989), estimated the expected bias 
of the observed objects induced by gravity. Their study, 
however, lacked definitive quantitative prediction due to 
the absence of a complete and reliable description of the 
matter correlation properties in the nonlinear regime. In 
this paper I would like to reconsider some of these calcu- 
lations but with precise quantitative results obtained in 
the quasi-linear regime. This is certainly not a complete 
derivation of the gravity induced peak correlations in a 
cosmological density field. In particular the positions of 
the cells are not constrained to be local maxima of the 



smoothed density field. The aim of my study is simply to 
derive the joint density PDF at two different distant loca- 
tions. In practice these results are applicable only to the 
largest mass concentrations observed in the Universe (the 
smoothing length should be larger than the correlation 
length) and for large separation. But it can give insights 
in the gravity induced large-scale biases. 

The paper is divided as followed. In §2 I present the 
general properties of the cumulants for joint densities in 
the quasi-linear regime and in the large-separation approx- 
imation. The quantitative derivations are given in §3. In 
§4 I discuss the validity of the large-separation approxi- 
mation. Eventually, §5 is devoted to comparisons with nu- 
merical simulations and §6 to a discussion of the possible 
astrophysical implications. 

2. General properties of the cumulants of the 
two-point density PDF 

I assume here that the cosmological density field evolved 
from Gaussian initial fluctuations following the gravita- 
tional dynamics of a pressure-less fluid in an expand- 
ing universe. Moreover the cosmological over-density field 
8(x)(= p(x)~p — 1) is assumed to be filtered by a given win- 
dow function so that the local rms fluctuations are small. 
Although the general properties discussed in this section 
are valid for any window function, most of the analytical 
results presented afterwards are valid for a top-hat win- 
dow function only. The numerical application will be done 
for an Einstein-de Sitter Universe. 

I assume here that the reader knows how cumulants, 
noted in the following (.) , are defined, otherwise see for 
instance Balian & Schaeffer (1989). 

2.1. Notations and diagrammatic representations 

The Fourier transform, 6(k), of the initial density field is 
defined by 

6{x) = J d 3 fc<S(fc)exp(ifc -x), (2) 

where 6(k) are random Gaussian variables. The power 
spectrum is given by 

(6(k)6(k')) = 6 Bliac (k + k') P(k). (3) 

The hypothesis of Gaussian variables is obviously of cru- 
cial interest for the derivation of the subsequent statistical 
properties of the density field. It implies in particular that 
the ensemble average of a product of an odd number of 
variables 6(k) is zero and that the ensemble average of a 
product of an even number is obtained by summing all 
the possible different products of cross-correlation when 
the 6(k) factors are associated by pairs. 

In the quasi-linear regime the evolution of the den- 
sity field is governed by the continuity, Euler and Poisson 



F. Bernardeau: The large-scale gravitational bias from the quasi-linear regime 



3 



equations. Let me consider and expansion of the grow- 
ing modes of the density field with the initial fluctuations 
field, 



8(x) = 8 (1 ^x) + 6 (2 \x) + 8 ( V(x) + 



(4) 



where 6^ is proportional to the initial fluctuations, 
quadratic, etc.. Following Goroff et al. (1986) one can 
write the expression of the density field at the order p, 
smoothed at the scale R, 



$(p) 



: J d 3 fci ...J d 3 fc p (5(fci) . ..6(k p ) 
x N v {k u ...,k v ) tF[|fc 1 + ...fc p | R]D p (t) 



(5) 




where D(i) is the time dependence of the linear growing 
mode, N p are homogeneous function of the wave vectors, 
W is the Fourier transform of the window function ap- 
plied to the field. (Note that this form is valid only for an 
Einstein-de Sitter Universe but admits simple generaliza- 
tions). 



6(p>= 



p lines 

Fig. 1. Diagrammatic representation of the term 6^{x). Each 
dot stands for a factor 8(k). 

In Fig. 1 I present a diagrammatic representation of 
the term h^ p \x). Each point appearing in the diagram 
stands for a factor 6(k); the lines for exp(ifc • x) and the 
hypothesis of Gaussian hypothesis implies that the en- 
semble average of an arbitrary number of such quantities 
should be calculated by associating the points by pairs (as 
in Figs. 2 and 3). 




Fig. 2. Example of a term contributing to (8 2 (xi)S 3 (a^)) • 
The circles at the left side of the figure are at the position X\ 
and the ones at the right side are at the position X2 ■ 



In the following I am interested in the calculation 
of the dominant order of the cumulants of the kind 
(b v (xi)b q (a^)) , where 6(x) is the smoothed local density 
at the position x. Such cumulants at their leading order 
in the quasi-linear regime are given by an expression of 
the form 



(6P(x 1 )6"(x 2 )) c 



decompositions i—1 i = l 



(6) 



where the decompositions are made in such a way that 
the integers pi and qj verify 



J2 Pi + J2 qi = 2(p+?)-2. 



(7) 



This constraint ensures that the number of lines connect- 
ing the factors is large enough for the term to be connected 
(terms in which the sum (7) is less than 2(p + q) — 2 do 
not contribute to the cumulants), but small enough that 
there are no extra lines that would lead to a negligible 
term in the quasi-linear regime. In Fig. 2 I present the 
diagrammatic representation of such a term contributing 
to the expression of ^<5 2 (a3i)<5 3 (a?2)) c - This particular term 
enters in the determination of 

(«( 1 )(i 1 )«( 3 )(zi) [^ 1 )(* 2 )] 2 ^ 2 )(a!2)) c . It is of interest 
to note that for q = (or p = 0) the constraint (7) reduces 
to 



E 



Pi 



2p-2, 



(8) 



and it corresponds to the calculations of the cumulants of 
the one-point density PDF. The latter are then expected 
to behave like (Fry 1984, Goroff et al. 1986, Bernardeau 
1992) 



s p (6\x)y- 



(9) 



where the power p — 1 is a direct consequence of the con- 
straint (8). The coefficients 5^ have been calculated ex- 
actly for a top-hat window function (Bernardeau 1994c). 
The calculation presented in the following is a generaliza- 
tion of these results to joint cumulants. 



2.2. The large separation approximation, factorization 
properties 

The general calculation of diagrams as in Figs. 2, 3 
is unfortunately quite complicated and, as far as I know, 
cannot be done in general. However, the calculations can 
be carried on in a limiting case, when the separation be- 
tween the cells is large. This limit can be understood in a 
diagrammatic point of view. The length of the lines joining 
two points in Fig. 2 is roughly proportional to the value 
of the two-point correlation function at the same scale. It 
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cell 1 



cell 2 



Fig. 3. Example of a term contributing to (8 2 (xi)S 3 (a^)) • 
The circles at the left side of the figure are at the position x\ 
and the ones at the right side are at the position X2 ■ 



ko only through the combination koR. In the large sep- 
aration limit the contributing value of ko in (12) will be 
of the order of 1/1*1 — a?2 1 , so that the function G(ko,R) 
has to be evaluated only in the limit ko — > 0. This makes 
the calculation quite simpler, since the cumulants (6) now 
read, 

(6P(x 1 )6^x 2 )) c = G p (0,R) G q (0,R) 

x J d 3 k P(ko) exp [ife • ( Xl - x 2 )] (13) 
= G p (0,R) G q (0,R){8( Xl )8(x 2 )). 
I then define the coefficients C PA by 

(6>>(x 1 )6«(x 2 )) e = C Ptq (S 2 ) p+q - 2 (6(x 1 )6(x 2 )) . (14) 



implies that long lines (for hierarchical, CDM like, power 
spectrum) give a negligible contribution compared to short 
lines. As a result, terms in which the number of long lines 
(the ones that connect the cells together) is minimal dom- 
inate the expression of the joint cumulants. For instance 
in Fig. 3 I present a diagram whose contribution to the ex- 
pression of ^<5 2 (a?i)(5 3 (£C2)} c is dominant compared to the 
one presented in Fig. 2 in the large separation limit. For 
two cells the number of lines joining the two cells should 
simply be one. The dominant contributions are now given 
by an expression similar to (6) but with 



"^2 Pi = 2p-l, 



(10a) 



It is then clear that the dominant term in the quasi-linear 
regime for the expression of C PA is constant in the limit 
of large separation (i.e., independent of d). Moreover from 
the previous equation the coefficient C p q is simply given 

by, 



_ G p (0,R) G q (0,R) _ ri 
/p\p /py 



(15) 



and is thus factorizable in p and q. This is one of the 
key properties obtained from the large-separation limit 
approximation. 



2.3. Examples of coefficients: C 2 \ and C 2 2 



2q - 1. 



(106) 



In this section I derive the values of C 2t \ and C 2<2 in the 
large-scale limit for a top-hat window function. In such a 
case the window function in Fourier space is given by 



So, in the large separation limit, only a subset of diagrams 
have to be calculated. The general form of such diagrams 
is given by 



6 p (x 1 )6 c '(x 2 )) c = / d 3 k G p (k ,R)G q (k ,R)P(k ) 



(11) 



exp [ifco • (xi — x j 



with 



G p (ko,R) = j d 3 fc! ...d 3 fc p _! 



x e n^.( fc *----< fc ^) ( i2 ) 

permutations {s} i = l 

x W [\k tl + ... + k Spt | R] Pih) . . . P(A p _i) 

with G\{ko,R) = 1 and where k Sj is either k\, fc p _i or 
ko- The determination of the integrals (11) is then reduced 
to the derivation of a single index series of integrals. The 
functions G p {ko, R), for dimensional reasons, depend on 



W(k) = — [sm(k) - k cos(fc)] 



(16) 



Before doing the calculations I recall few properties of 
this window function, given by Bernardeau (1994a), 



1 

47T 



and 



dfi W[\ki + k 2 \] 



-W(h)W(h), 



(ki-k 2 f 



h 2 h 2 



— / dnw^Ofei + fe 2 |i i l 



k\ ■ k 2 



W(k 2 ) + -k 2 W'(k 2 ) 

o 



(17) 



(18) 



where the integration is made indifferently over the angu- 
lar part of k\ or k 2 . 
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C'2,1 



The expression of the coefficient 6*2,1 is given by 
4 



Using again the properties (17, 18) we get 



: J d 3 fc d 3 fci P(k ) P{ki) exp [ifc • (*i - x 2 
:W(k R) Wilko + kilRlWihR) # 2 (fc ,fci) 



(19) 



{8( Xl )8( X2 )){8i) 

: J d 3 fc d 3 fci P(k ) P(h) 

: exp [i(fc + fci) • (xi - x 2 )] 

:W(k R) Wilko + kilRlWihR) F 2 (fc ,fci) 



with 



5 1 fc • fci 1 fc • fci , 2 fk ■ ki 
N 2 (k , fci) = - + 77^2— + 7j— p- 



7 2 H 



7 V fco fci 



The second term appearing in (19) is in fact negligible 
since it involves the presence of two lines joining the cells 
at the positions x\ and x 2 . The first term of this expres- 
sion can be calculated using the properties (17-18) (inte- 
grating first on the angular part of the wave vector fci). 
It leads to the expression, 

_ 68 , l dlog(£ 2 ) l dlog(%i)% 2 )) 
62,1 "21 + 3 dlogi? + 3 dl^R • (20) 

In the large separation limit the logarithmic derivative 
of (6{xi)6(x 2 )) with respect to the smoothing scale is ex- 
pected to be zero: this cross-correlation is simply the value 
of the two-point correlation function of the un-filtered field 
for the separation \x\ — x 2 \. The resulting value for C 2t \ 
in the large separation limit is simply given by* 

_68 l dlog(^) 

21 + 3 dlogi? • (21) 

Note that this result depends on the local index of 
the two-point correlation function only at the smoothing 
scale, not at the separation at which the cells are. But this 
is only true for the large separation limit. The corrective 
term appearing in (19) for instance is likely to have a more 
complicated dependence with the spectrum. 

The derivation of the coefficient C 2t2 can be done di- 
rectly from the expression (valid in the large separation 
limit), 

c 16 

(-^2,2 



<^<%i)% 2 )) 

x j d 3 fcid 3 fc d 3 fc 2 Pih) P(k ) P(k 2 ) 

x exp [ifc • {x 1 - x 2 )] W(k 1 R) W(k 2 R) 
x W[\ki +k \R] W[\k 2 + k \R] 
x F 2 {k 0) ki)F 2 {k 0) k 2 ) 



(22) 



* This result is in fact related to the value Q given by 
Fry (1984, eqn. 38) since one should have C 2) i = 2Q. 



C 2 



2 68 dlog(<S 2 ) 1 / dlog(<S 2 )' 



3 21 dlogi? 9 \ dlogi? / 

I 2T / d3fc2jP ( fc2 ) W ^ 2 R ) k ' 2 R W '^ 2 R ) 
x exp[ifc 2 • (xi - x 2 )]/(8(xi)8(x 2 )) 

d 3 k 2 P(k 2 ) W(k 2 R) k 2 R W'(k 2 R) 



ldlog<£ 2 ) ,,3, 



9 d log R 

x exp[ifc 2 • (xi - x 2 )]/(8(xi)8(x 2 )) 

+ i J d 3 k 2 P(k 2 )k 2 R 2 [W'{k 2 R)f 

x exp[ifc 2 • (xi - x 2 )]/(8(xi)8(x 2 )) 

The last three terms of this expression are negligible in 
the large separation limit due to the presence of the factor 
k 2 R which, for k 2 ps 1/1*1 — x i\, is small. The resulting 
value of C 2)2 is then simply given by 



C 2 



r 2 

°2,1- 



(23) 



It confirms for this particular example the factorization 
property (15). 

3. The generating function of the joint cumulants 

The quantity we are interested in is actually the generat- 
ing function (3{y\, y 2 ) defined by 



00 v q 



(24) 



p = l,q = l 

Due to the factorization property (15-16) the derivation 
of the C Pt q series can then be done with the derivation of 
C Pt \ coefficients only, so that 

h(yi,V2) = P(yi)P(y2) (25) 

with 



P(y) = E^.i 



P =i 



(26) 



3.1. When the smoothing effects are neglected 

When one neglects the effects of the window function that 
is by assuming that 

W[\k 1 + ... + k p \R}& W[klR]. . W[k p R] in (12) the 
derivation of these integrals is quite simple. The factors 
C Pt \ can be seen as sums of products of vertices v p , geo- 
metrical averages of the function N p . These products cor- 
respond to tree summations for trees having one "free leg" . 
Note that the S p parameters (Eq. [9]) on the other hand 
are obtained with tree summations for trees without any 
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free leg. These tree summations have been considered by 
Bernardeau & Schaeffer (1992) and Bernardeau (1992). 
They are all related to the generating function of the ver- 



tices v n 



In the following the scale R\ is identified to the scale R 
of the filtering scale and the scale R 2 is identified to the 
typical distance between cells, assumed to be larger than 
the smoothing scale. 



2^ v p 

p = 



(-r)P 



pi 



(27) 



From Bernardeau (1992) the function Q(t) has been shown 
to describe the spherical model, and is thus a known func- 
tion for any cosmological model. The generating function 
f3(y) is then given by 



(28) 



From this implicit relationship one can easily relate 6*2,1 to 
the value of v 2 = 34/21 (for an Einstein-de Sitter Universe, 
and one find, 



6*2,1 = 2^ 2 



68 
21' 



(29) 



which corresponds to the case (21) when n 



3.2. When the smoothing effects are taking into account 

Taking into account filtering effects is of crucial impor- 
tance for a practical point of view. Since the pionneer- 
ing work of Goroff et al. (1986) the filtering effects are 
known to change the values of the S p parameters. A sim- 
ilar change is going to affect the coefficients C VA . The 
results of the previous sections show that it is indeed the 
case for the 6*2,1 parameter. It is of interest to note that 
the relationship between 6* p ,i and S p that holds in the 
unsmoothed case does not hold anymore. For instance we 
have 6*2,1 = 2/3 S3 for the unsmoothed case but this is not 
true when the smoothing is taken into account (it would 
have given 68/21 - 2(n + 3)/3 instead of 68/21 - (n + 3)/3 
for 6*2,1). The proper calculation of the coefficients 6* p ,i 
should then be done with care. 

3.2.1. The case of two concentric cells 

The basic property I use to do this calculation is to re- 
mark that the parameters 6* p ,i are also involved in the 
expression of i^6 R (x)6r 2 (x)) , where the smoothed den- 
sity fields are taken at the same position but when the 
scale R 2 is very large compared to R\, R\ <C R 2 . This 
property is only true for the large separation limit. The 
reason is that for the calculation of such a quantity the 
factor exp[ifco • (*i — a^)] W[&o R] in (13) has simply been 
replaced by VF[&o R2], so that this new cumulant is given 

by 



$ p rM 



G p (0,i?i)(%^ 2 ) 
c pa( 6 r 1 ) P 1 {SrJr 2 



(30) 



3.2.2. The joint PDF in Lagrangian space 

The derivation of the 6* p ,i parameters follows then the 
same scheme than the one developed for the S p param- 
eters. The fundamental property is that the filtering in 
Lagrangian space (that is for a given mass scale) does not 
change the values of the cumulants. This property had 
been discussed for the case of a unique cell but can be 
easily extended to any number of concentric cells. The 
reason is that this property is due to angular properties 
of the top-hat window function similar to (17-18) (actu- 
ally to the properties (B3) and (B4) of Bernardeau 1994a) 
that involve only the angular part of the wave vectors. 
As a result the filtering effects simply factorized away for 
concentric cells, i.e., the effect of filtering appears only 
in the value of (6 R ) (which obvioulsy depends on R), of 
(6r 1 6r 2 )> but not in the values of the C VA coefficients. The 
derivation of the joint density PDF in different concentric 
cells can then be performed using the technics presented 
by Bernardeau & Schaeffer (1992), with the generating 
function of the vertices given by Q(t) in Eq. (27). 

So, let me now consider the joint density PDF in 
two different concentric cells of two different mass scales, 
M\ <C M 2 ■ Following the calculations made by Bernardeau 
k Schaeffer (1992) we get 

-If dt/i dt/2 



P L (61,62) 



with 



X L (yi,V2) 



4tt 2 J <T 2 (Mi) <t 2 (M 2 ) 

I" l / n , vi6i 

exp X (2/1,2/2) 



(31) 



2/2^2 



<r 2 (Mi) [ 



,7 2 (Ml) <7 2 (M 2 ) 



<r 2 (M 2 ) 



V2Q{t2) 

a 2 (M 1 ,M 2 ) 
<r 2 (M 2 ) 



1 L d r , L, 



(32) 



y2 Tr Q(rt) 



y2 TT Q(r!f) 



<r 2 (M lt M 2 ) d L 



where <r(Mi), <r(M2) and a{M\,M'f) are respectively the 
rms density fluctuations at the mass scale M\ , at the mass 
scale M2 and the cross correlation between the mass scales 
M\ and M 2 - The superscript L designs quantities related 
to Lagrangian space calculations. Note that the density 
PDF at the mass scale M\ is given by 



P L (6i) 



<iy 

2ma 2 (M 1 ) 



exp 



X L (2/i,0) 



2/i^i 



<r 2 (Mi) 



(33) 
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Our present objective is only to get the generating func- 
tion of the coefficients C Pt \. It is thus necessary to limit 
ourself to the case where only the cumulants (^r 1 ^R 2 ) c 
and (i5r 2 ) are not zero. In such a case the function 
X L (2/1,2/2) reads, 

a 2 (M u M 2 ) L 

(34) 



x L (1/1,1/2) =x L (i/i,o) + 



a 2 (M 1 )a 2 (M 2 ) 



2a 2 (M 2 ) 
where T\ is given by 



(35) 



In Lagrangian space the generating function f3(y) is simply 
given by 



p L (y) = r L (y). 



3.2.3. The joint PDF in Eulerian space 



The problem is then to relate the Lagrangian density PDF 
to the Eulerian density PDF. The first remark of interest 
is that at the mass scale of M 2 no differences are taken into 
account: the reason is that the cumulants we are interested 
in involve only the linear order at the mass scale M 2 ■ The 
Lagrangian mass scale M 2 and the Eulerian scale R 2 are 
then identical and the actual value of 8 2 can be seen as 
a simple external constraint. We are then confronted to 
a problem similar to the one considered by Bernardeau 
(1994c) and a similar development will be followed. The 
relationship between the density PDF in Eulerian scale 
(<5i , ^2) and the one for Lagrangian scale Mi is 
simply that (see Bernardeau 1994c) 



'So 

with 



d«i (1 + «i)p£ («i,« 2 ) 



3M 



d6 lP % Io (6 1 ,6 2 ) (37) 



(38) 



The derivation of the cumulants (in the large scale limit) 
can be done by identifying the exponential factor in the 
expression of p E (6i, 6 2 ). It is then interesting to note that 
integrating and differentiating are not going to affect the 
expression of the exponential term so that a simple iden- 
tification can be done, taking advantage of the mass-scale 
relationship (38). 

The calculation of the expression of p L (6i , 6 2 ) can be 
done by a saddle point approximation. The saddle point 
position is given by the system of equations, 



d 

dyi 
d 

dy 2 



X L (2/1, 1/2) 



x L (2/1, 1/2) 



<5i 2/1 



<t 2 (M0 

62 1/2 



<7 2 (M 2 ) 



(39) 



The solution, (2/1,2/2), of this system is given by, 
Then, at the saddle point position we have, 



(40) 



x L (y{,y 2 ) + 



1/^1 



y 2 h 



<7 2 (Mi) <7 2 (M 2 ) 

if(y°) <r 2 {M u M 2 ) 



2<r 2 (M 1 ) <7 2 (Mi)(7 2 (M 2 ) 



*2 



(41) 



2<r 2 (M 2 )' 



when the expression is written up to the linear order in 
a 2 (Mi, M 2 ). The value of y® is the one given by the saddle 
point approximation for the calculation of p L (Mi) in (33), 
so that, 



(36) g[r 1 L (y° 1 )] = 8 1 . 



(42) 



In Eulerian space a similar expression is expected so 



that 

x E (yl,y 2 ) + 



y 2 h 



a 2 (R x ) a 2 (R 2 ) 
rffo?) a 2 (R u R 2 ) 



(43) 



2a 2 (R 1 ) a 2 (Ri)a 2 (R 2 ) 

with 

G E [r 1 E (y° 1 )] = h. 



62p E (y° 1 ) 



2a 2 (R 2 ) 



(44) 



The expressions of Q E , t e and j3 E can then simply be 
obtained by identification. For that the mass scales M\ 
and M 2 have to be transformed into the physical scales 
Ri and R 2 . For the mass scale M 2 , it is rather easy since 
this scale is supposed to remain linear. As a result we 
simply have 



M, 



47T o 



(45) 



The mass scale M\ and the scale R\ are related to each 
other with a relation similar to (38), 

47T 



M x =7,(1 + 6!) ^-Rf. 
The identification of the various terms then gives 



(46) 



r E (yl) = r L (yl) 



g E (r E ) = g(r L ) = g 

f3 E (y 1 ) = h(y° 1 ) 



a(R 1 [l + g L (T E )y/ 3 ) , 
^(R^i + g^)]^) 



t(Ri) 



(47) 



t(Ri) 



T(i?i[l + a i (r i )] 1 /3) 



For this result we can neglect the possible variations of 
<t(R\,R 2 ) with that is 



<t(Ri[1 + G E (t e )] 1 I\R 2 ) ^a(R u R 2 ) 



(48) 
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since R\ <C R 2 - This is a simple generalisation of the re- 
mark leading to (21) from (20) for the calculation of the 
coefficient 6*2,1. 

It is interesting to write the results in terms of the 
function T E (y), 



P E (v) 



a(R) 



<r(i?[l + a £ (r £ )] 1 /3) 



(49) 



The function j3 E (y) is the generating function of the coef- 
ficients C Pt \ as defined in Eq. (26). Although it has been 
obtained from the statistics of two concentric cells, it is 
identical to the one corresponding to two distant cells of 
the same radius. Note that, unlike the Lagrangian case, 
the function j3 E is no more equal to the function t e as 
a naive extrapolation from the unsmoothed case would 
have given. The relation (49) is the central analytical re- 
sult of this paper, and all the applications that followed 
are derived from this formula. 



3.2.4. Examples of coefficients C 1 



To start with it is possible to get, with the expansion of 
j3 E (y) with respect to y, the first coefficients C Pt i. As an 
illustration I give here the first six, 



C 2,i-2i+ T , 



(50) 



C31 



11710 61 



441 



2 _l 72 



C 4 ,i =353.1 + 205. 171 + 38.67 7l 2 + 2.33 7 f + 12.03t 2 

+ 2.22 7l72 - 0.22 T3 ; 
C* 5 ,i =6511.17 + 5497. 171 + 1715. 7 7l 2 + 233. 5 7 f + 11. 67? 

+ 390.2772 + I6O.277172 + 16.1l7 2 7 2 + 1.852 72 

+ II.5773 + 2.3467173 + O.I23474; 
C* 6 ,i =153708 + 17105171 + 754607 2 + 164607f 

+ 1770 7 f + 74.7 7l 5 + 1355172 + 88787172 

+ 1916 7l 2 72 + 135.7 7 f72 + 172.7 7 2 + 36.9 7 i7 2 2 

+ 502.873 + 217.4 7 i73 + 23.21 7l 2 73 + 4.2O7273 

+ 8.9474 + I.97174 + O.O675; 



with 



1p 



dP\og(6 2 ) 
dlog p R ' 



(51) 



The expression of 6*2,1 have already been obtained from 
direct calculation in (21). I also checked that the coef- 
ficient 6*3,1 can also be re-derived directly from formula 
(12) and using the properties (17, 18) of the top-hat win- 
dow function. The results are given here for an SI = 1, 
A = Universe, but, as for the S p coefficients (Bouchet et 
al. 1993, Bernardeau 1994b), the £1 and A dependence of 
these coefficients are expected to be small. 



3.3. The bias function 

The interest of such results is that it is then possible to 
get the expression of the joint density PDF in two distant 
cells. In the following I assume that the two cells have the 
same radius R. The rms density fluctuations in such cells 
will be denoted <r{R), and the cross-correlation between 
the two cells, a 2 (d). Moreover I will suppress the E super- 
script but it is always assumed. The expression of such a 
joint distribution function can be obtained from the form 
(31), with the function exp[x(j/i, J/2)] in (32) calculated in 
the large separation limit that is up to the linear order 
in a 2 (d). In such a limit it reads (Bernardeau & Schaeffer 
1992) 



exp[x(2/i,2/ 2 )] 



1 



exp 



aHR) 



^(2/2) 



<T 2 {R) <7 2 (i?)_ 

As each term appearing in (52) can be factorised in 

p{8 1 ,8 2 ) = p(6 1 ) P (6 2 ) + a 2 (d) p 6 («i)p 6 (« 2 ) 

with 

Ay -P(y) 



(52) 



(53) 



Pb(S) 



27ri a 2 (R) a 2 (R) 

f(y) , 6 y 



exp 



(54) 



a 2 (R) a 2 (R)\ 
It is then quite convenient to define the bias function 



b(8) 



Pi(6) 



/dy/?(y)exp + 



<p(y) 



6 y 



P( 6 ) J dy exp 

which implies that the joint density PDF reads, 



(55) 



p(6 1 ,6 2 )=p(6 1 )p(6 2 ) [l + <7 2 (i?i,i? 2 )6(«i)6(« 2 )] . (56) 



3.3.1. Normalisation properties 

In the following I will discuss the various properties of the 
function b(8). The first remark is that 



d6 b(6)p(6) = 



d6 8 b(8) p(8) = 1. 



(57) 



(58) 



In Fig. 5 I present the shape of the bias function for a 
particular set of parameters. The general features that are 
seen is that it is a growing function of 8. It is (roughly) 
of the sign of 8. Note that for a purely Gaussian field, one 
would have 



^Gaussian 



(59) 
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Fig. 4. The approximate form (64) (thin solid line) compared 
to the exact integration of (55) (thick solid line). The thin 
dashed line is the form (65) and the thick dotted line is the 
Gaussian case (59). 



This 



is shown as a thick dotted line. 



3.3.2. Analytic approximate forms 

It is useful to remember that the function Q{t) in (47,49) 
can be fruitfully approximated by 



2t 



-3/2 



(60) 



In the following I will also assume that the power spectrum 
can be approximated by a power law behaviour so that 



P(k) oc k n . 
It implies that 

and 



(61) 



(62) 



pE = - (1 + 0£)(n+3)/3 | (1 + gEy2/3 _ ^ (gg) 

An analytical approximate formula for b(8) can be ob- 
tained with the saddle point approximation. It implies 
that 



b(6) 



" 3 (l + 5 )(«+3)/3 [ (1+(5) -2/3_ x 



(64) 



In fig. 4 I show this approximate form compared to an ex- 
act integration of the expression (55). the approximation 
is very good for 6 £ but becomes quite inaccurate when 
the density contrast is large. 



Another possible approximation is the following ap- 
proximation, 



b(8) 



Ay [] E (y)exp[(l + 6)/* 2 } 
I [ dy p(j/)exp[(l + «)/<7 2 ], 



(65) 



which is a function of (1 + 6) /a 2 only. As shown in Figs. 
4 and 7, this approximation is valid when the value of 6 is 
large, and the larger a is the better it is. 




Fig. 5. The variation of the bias function for different values 
of a. The approximate form (64) (thin lines) are compared to 
the exact integrations (55) (thick lines). The dashed line, solid 
line and dot dashed line are respectively for a = 0.3, a = 0.5 
and a = 0.8. 




1 + 5 

Fig. 6. The variation of the bias function for different values 
of the power law index n. The solid line, dashed line and dot 
dashed line are respectively for n = — f, n = — 2, n = — 3 
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Fig. 7. The bias function expressed as a function of (1 + 8) /a 2 
for different values of a. The solid line is the approximation 
(65). The dashed line, long dashed line and dot dashed line are 
respectively for a = 0.3, a = 0.5 and a = 0.8 and the thin lines 
correspond to negative values of b(8). 
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Fig. 8. The cumulative bias function (66) expressed as a func- 
tion of (1 + 8) I a 2 for different values of a. The symbols are the 
same than in Fig. 7. 



In Figs. 5 and 6 I present the shape of the function b(6) 
for different values of the index n and of a. the variation 
of the shape of b(8) for different values of n are impor- 
tant for negative values of 6. For positive values of 6 the 
variations are not so important. The approximation (64) 
is however quite inadequate to see this latter property. 
The approximation (65) would have been more accurate 
for this regime. 

The Fig. 7 is quite instructive to show that the a vari- 
ations, for large values of 6, are well described by the ap- 
proximation (65). The function b(8) indeed appears to be 
dependent only on (1 + 6) /a 2 . 



A quantity that I will consider for the comparison with 
numerical simulation is the bias function for a threshold. 
It is defined by 



£ as b(6) P (6) 



(66) 



In Fig. 8 I give the shape of the function 6(> 6) expressed 
as a function of (1 + S)/u 2 and for different values of u. 

3.3.3. Asymptotic behaviour 

When 6 is large the form (54,55) naturally reduces to the 
form (65). The integrals when 6 is large are dominated by 
the singular point y s of the expression of <p(y), 



(p(y) - (p s « -a s (y - y s ) 3/2 . 

At this point the function j3 E (y) reads 

p E (y)-pf *-b s (y-ys) 1 ' 2 . 



(67) 



(68) 



The asymptotic behaviour can then be calculated follow- 
ing the method of Bernardeau & Schaeffer (1992) and 
Bernardeau (1994c). We obtain 



b{6) pa b* 



with 



3 6, 
2 a. 



(69) 



(70) 



Peculiar values of these parameters are given in table 1 
for different values of n . 

Table 1. The parameters of the critical point as a function 
of the spectral index, n, for the bias function. 



n 


a s 


b s 


&* 


-3 


1.84 


0.528 


0.431 


-2.5 


2.21 


0.644 


0.438 


-2 


2.81 


0.836 


0.447 


-1.5 


3.93 


1.206 


0.460 


-1 


6.68 


2.154 


0.484 



The resulting value of is slightly lower than 0.5 and 
is only weakly dependent of n. It is also interesting to 
note that this asymptotic behaviour gives a lower value of 
the bias than for the Gaussian case for a given value of 
6 /a 2 . Note, however, that in the present paper the value 
of 6 is not the initial overdensity of the fluctuation but its 
final value. The two appraoches, peak correlation in the 
initial Gaussian density field and correlation in the final 
quasi-linear field, then cannot be compared directly. 
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4. Validity of the large separation approximation 

The validity of the large separation limit is quite difficult 
to test in general. The reason is that a lot of terms have 
been neglected to get the results (50) for the large sep- 
aration limit. With numerical integrations it is however 
possible to estimate the value of 6*2,1 for any separation. 
To do the numerical calculations I use the shape of the 3- 
point correlation function obtained from the second-order 
perturbation theory given by Fry (1984), 



10 

y 



4 -3 - 2n + (n + 3) 2 cos 2 
7 n 2 
+ cyc.(l,2,3) 



3 + n ( r 21 r 3i 

cos(0i) 1 

n \ r 3i r-2i 



6(^1,^2) h(ri,r 3 ) 

(71) 



where 9\ is the angle between the vectors r 2 — r\ and 
r 3 — r\ and r 2 \ and r 3 \ are the lengths of these vectors. 
Two points are then assumed to be in a cell of radius R 
at the position xo and the other in another cell of radius 
R at the position x\ with |a?o — x i\ = d. The resulting 
coefficient, 6*2,1, is then a function of d/R and n only. In 
Fig. 9, I give the results for power law spectra of index 
n = —1.3 and n = —1. 



3.5 



2.5 




d/R 



Fig. 9. Numerical evaluation of the 6*2,1 coefficient as a func- 
tion of separation d. The calculations have been made for 
n = —1.3 (thin lines) and n = — 1 (thick lines). The circle 
is the d = prediction (in such a case 6*2,1 (d) = S3) and the 
dashed line is the large separation limit (21). 



When d = the two cells are exactly superposed 
so that the coefficient C*2,i is exactly S3 (Eq. [9]). For 
d/R 1 one expects the value of C*2,i to converge to the 
large separation limit (21). And it can be seen that this 
limit is indeed obtained very rapidly, as soon as d/R 3. 
Actually the various corrective terms (appearing in (19) 
and (20), tend to cancel each other. It is a bit surprising 



however that C*2,i is not a monotonic function of the dis- 
tance. It has a minimum for d ps 2R, that is when the 
cells are side by side. Actually although the values of S3 
and 6*2,1 in the large separation limit depend only on the 
local index at the smoothing scale the whole shape the 
function C*2,i as a function of the distance may depend on 
other characteristics of the power spectrum. 



5. Comparison with numerical simulations 

For comparison with numerical results, I use the nu- 
merical simulation done by Couchman (1991) with CDM 
initial conditions for h = 0.5. This simulation uses an 
adaptive P 3 M code with 2.1 10 6 particles. The simulation 
was made in a cubic box of 200/i -1 Mpc size with periodic 
boundary conditions. I used the simulation at a time-step 
for which the linearly extrapolated rms density fluctuation 
in a spherical box of 8/i -1 Mpc radius is unity. 



5.1. The measure of the coefficient C 2,1 

For the determination of the various statistical properties, 
I measured the number of particles in 200 3 spherical cells 
disposed on a regular grid. Each sphere is then associated 
by pairs to 6 neighboring spheres that are at the distance 
d along each direction of the grid. By summing over the 
particle content of each sphere (appropriately weighted) 
multiplied by the content of its neighboring spheres, it 
is possible to estimate the value of statistical quantities 
such as the coefficients C VA . However, an accurate mea- 
surement of merely the coefficient C*2,i turns out to be 
quite difficult to obtain even for rather small separation. 
For instance I found 



C*2 



1.9 ± 1.2 



(72) 



for R = 10/i -1 Mpc and d = 20/i -1 Mpc where the error- 
bar has been obtained by dividing the simulation in four 
subsamples. As it can be noticed the estimated errors are 
pretty large, making the comparison with the theoretical 
prediction rather un-conclusive. A precise analysis of this 
problem is beyond the scope of this paper. Note that the 
possible sources of errors for the one-point density PDF 
have been investigated in details by Colombi et al. (1995). 
It turns out the such coefficients are very sensitive to finite 
size effects, a problem which might be amplified here by 
the fact that I considered pair of cells. It may also be 
due to the fact that the number of pairs that have been 
used for the determination of joint moments is slightly too 
small, making the sampling of the large density tail of the 
distribution too poor for an accurate determination of the 
moments. These problems led me to focus my analysis on 
the behaviour of the global bias function. 
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0.1 1 10 

1 + 5 




1 + 6 

Fig. 10. The measured bias function for R = 5/t -1 Mpc in a numerical simulation with CDM initial conditions. For this scale 
a = 0.95 and n = —1.3. The upper panel shows the measured density distribution p(8) compared to the predictions from 
Perturbation Theory (solid line). The central panel shows the measured function b(> 8) with the formula (73) (squares) and 
with the formula (77) (triangles). The distance between the centers of the cells is d = 10/t -1 Mpc. The bottom panel shows the 
function b(> 8) obtained with the formula (77) when the distance between the centers of the cells varies from d = 10/t _1 Mpc 
to d = 20/t -1 Mpc. The error-bars have been obtained by the measure of these quantities in four different subparts of the 
simulation. The solid lines are the prediction of Perturbation Theory for the large separation limit, Eq. (55). 
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Fig. 11. The measured bias function for R = lOh 1 Mpc. For this scale a = 0.46 and n = — 1 and the distance between the cells 
is d = 20/t -1 Mpc. The symbols are the same as in Fig. (10) 



5.2. The measure of the bias function 

The direct measure of the bias function turns out to be 
easier to do in practice. I considered two ways to get the 
bias function. In the first method I estimated 



d«i / d^ 2 (1 + <5 2 )p(<5 1 ,<5 2 ) 



1 



E 



200 3 ^ 6 

cells N(i)>N 



(75) 



i 



(6162 



f~d6 1 f~d6 2 (l + 6 2 )p(6 1 ,6 2 ) 
P(> 6 ) 



(73) 



d6 1 / d6 2 (l + 6 1 )(l + 6 2 )p(6 1 ,6 2 ) 



200 



N 2 (i) 
6 ; 



(76) 



where 
(61S2) 



d8 1 d8 2 (1 + ^)(1 + 6 2 )p(6 1 ,6 2 ) - 1. (74) 



where N(i) is the number of particles in cell i and N 2 (i) is 
the total number of particles in the six neighboring cells 
at the distance d. 

The second method is based on the determination of, 



b 2 (> 8 ) 



1 



(6162 



J?° d<5i d6 2 p(6 1 ,6 2 ) 



1 



1/2 



(77) 



In practice the quantity p(> 60) is estimated by deter- 
mining the fraction of cells in which the number of par- 
ticles is greater than No = (1 + 8o)N (N is the mean This expression is obtained by counting the fraction of 
number of particles in a cell). Then the integrals interven- pairs of cells for which the number of particles exceeds a 



ing in the expression of &i are estimated by 



fixed number Nn in both cells. 
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The main interest in considering these two approaches 
is that for the form (56) the two functions should coincide. 
It is thus a test for the factorization property (56). 

The numerical results are given in Figs. 10 and 11. 
The comparison of the bias functions (73) and (77) (re- 
spectively the squares and the triangles) show that they 
coincide nicely (within the error-bars). This is an impor- 
tant result since it means that the form (56) for the joint 
density PDF is very robust and holds even when the large 
separation limit is not reached. It is of interest however to 
have in mind that the properties (15) and (56) are equiv- 
alent only in the large separation limit. So this numerical 
result does not mean at all that the property (15) holds 
for any separation. 

The dependence of the bias function with the sepa- 
ration is presented in Fig. 10. For the distance at which 
the bias function can be reliably measured it is found to 
be very stable and independent of the distance. It checks 
another aspect of the form (56), that is that the depen- 
dence with the distance d enters only through the value of 
the mean two point correlation function between the two 
cells. It implies in other words that the bias is very weakly 
dependent of the distance. 

The Figs. 10 and 11 also show that the theoretical pre- 
diction for the bias function gives a very good description 
of the measured shape of 6. This is quite remarkable since 
the measures are not made at a very large separation. 

6. Discussion 

A reliable description of joint matter density PDF-s has 
been obtained based on PT. The analytical results are 
valid in principle in the large separation limit and for a 
smoothing scale larger than the correlation length. Com- 
parison with results of numerical simulations proved that 
they provide us with a good description of the joint density 
PDF for a broad range of parameter. In particular they are 
still accurate when the rms density fluctuation approaches 
unity at the smoothing scale and when the distance equals 
only two times the distance of the smoothing length. 

Remarkable properties have been obtained using PT in 
the large separation limit. In particular the property (15) 
of the C Pt q coefficients leads to the factorization form (56). 
The latter property has been confirmed by the numerical 
analysis and is shown to be in fact a good description of 
the joint density PDF for any value of the distance (larger 
than 2R). This implies that the dependence with the dis- 
tance enters only through the cross-correlation between 
the cells, u 2 (d), and that the dependence of the bias factor 
with the two densities can be factorized in 6(61)6(62). It 
is however worth reminding that the extension of the fac- 
torization property (56) to small distances does not prove 
that the property (15) also holds for small distances. 

Moreover, the analytic work gives the expression of the 
bias function b(8) from the large separation limit approx- 



imation, and these analytical results are in good agree- 
ment with the numerical measurements. It suggests that 
the gravitationally induced bias does not vanish at large 
scale and is well described by the quasi-linear dynamics. 

The accuracy of the theoretical prediction is very en- 
couraging and shows that the theoretical expression of the 
bias can indeed be used for quantitative predictions. 

Its implication for the estimation of the error-bars on 
statistical measurements related to the one-point density 
PDF due to finite sample effects will be exploited in a 
coming paper. 

For direct astrophysical applications, (e.g. the deter- 
mination of what could be the bias associated with the 
observed clusters) the analytic results obtained in this pa- 
per have two important limitations; 

• contrary to the calculations of Bardeen et al. (1986) 
this is not the correlation of peaks. However it is ex- 
pected that, in the rare event limit, the results are not 
affected: a very dense spot in the density field is indeed 
expected to be a peak. 

• The other limitation is that the calculations have been 
made for a large smoothing length, larger than the size 
of the clusters. This property is due to the use of PT 
calculations and can hardly be overcome except with 
pure numerical analysis. 

It anyway indicates that gravity is expected to induce 
specific biases, that are not only due to peak selection 
effects as for the results of Bardeen et al. (1986). Note 
that we anyway do not expect qualitative changes from 
this previous picture: the bias factor associated with the 
clusters is expected to be independent of scale and propor- 
tional to its mass in the large richness limit (eq. 69). Such 
properties were also predicted to be a generic result for 
the correlation of dense spots in a strongly nonlinear cos- 
mological field by Bernardeau & Schaeffer (1992). I leave 
for another paper the estimation of more quantitative pre- 
dictions. 
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